#include <iostream>
#include <math.h>

int sign(double x)      //符号函数
{
    if (x > 0)
        return 1;
    else if (x == 0)
        return 0;
    else
        return -1;
}

void swap(double &x, double &y) //交换函数
{
    double z = x;
    x = y;
    y = z;
}

// 目标函数的抽象类
class Function
{
    //diff可采用差商代替导数
public:
    virtual double operator()(double _x) = 0;       //求函数值
    //virtual double derivetive (double _x) = 0;      //求导数值,不是必需
    // = const{return ((*this)(x+_epsL)-(*this)(x-_epsL)) / (2*_epsL)}
};

//用户输入函数
class MyFunc : Function 
{
private:
    double (*fp)(double);
    double (*dfp)(double);
public:
    //fp为用户输入函数f(x)的指针,dfp为导数
    MyFunc(double (*_fp)(double)){fp = _fp;}    //允许只输入函数
    MyFunc(double (*_fp)(double), double (*_dfp)(double)){fp = _fp; dfp = _dfp;}    //输入函数和导数
    double operator()(double _x){return fp(_x);}
    double derivetive (double _x){return dfp(_x);}
};

// 求解器的抽象类
class EquationSolver
{
public:
    virtual double solve() = 0; //纯虚函数是通过在声明中使用 "= 0" 来指定
};

class Bisection : public EquationSolver
{
private:            //可以考虑打包成structure
    double a;       //left limit
    double b;       //right limit
    double delta;   //distance of x
    double eps;     //error of f(x)
    int M;          //max step
    MyFunc f;      //MyFunc为用户定义
public:
    Bisection(double _a, double _b, double _delta, double _eps, int _M, MyFunc _f);
    ~Bisection(){};
    double solve();
    void condition();       //show elements in class
};

//Bisection's member function
Bisection::Bisection(double _a, double _b, double _delta, double _eps, int _M, MyFunc _f):
a(_a), b(_b), delta(_delta), eps(_eps), M(_M), f(_f){}

double Bisection::solve()
{
    double h = b - a;
    double u,c,w;
    for(int k = 1; k < M; k++){
        u = f(a);
        h /= 2;
        c = a + h;
        w = f(c);
        if(fabs(h) < delta || fabs(w) < eps)
            break;
        else if (sign(w) == sign(u)){
            a = c;
        }
    }
    return c;
}

void Bisection::condition()
{
    std::cout << "a = " << a << "\t" << "b = " << b << "\t" << "M = " << M << std::endl;
    std::cout << "delta = " << delta << "\t" << "eps = " << eps << std::endl;

}

class Newton : public EquationSolver
{
private:
    double x;       //x0
    double eps;     //error of f(x)
    int M;          //max step
    MyFunc f;      //MyFunc为用户定义
public:
    Newton(double _x, double _eps, int _M, MyFunc _f);
    double solve();
};

Newton::Newton(double _x, double _eps, int _M, MyFunc _f):
x(_x), eps(_eps), M(_M), f(_f) {}

double Newton::solve()
{
    double u;
    int k;
    for(k = 0; k < M; k++){
        u = f(x);
        if(fabs(u) < eps)
            break;
        x -= u/f.derivetive(x);
    }
    //std::cout << "k = " << k << std::endl;
    return x;
}

class Secant : public EquationSolver
{
private:
    double x0;       //x0
    double x1;
    double delta;
    double eps;     //error of f(x)
    int M;          //max step
    MyFunc f;      //MyFunc为用户定义
public:
    Secant(double _x0, double _x1, double _delta, double _eps, int _M, MyFunc _f);
    double solve();
};

Secant::Secant(double _x0, double _x1, double _delta, double _eps, int _M, MyFunc _f):
x0(_x0), x1(_x1), delta(_delta), eps(_eps), M(_M), f(_f){}

double Secant::solve(){
    double u = f(x1);
    double v = f(x0);
    double s;
    int k;
    for (k = 2; k < M; k++){
        if (fabs(u) > fabs(v)){
            swap(x0,x1);
            swap(u,v);
        }
        s = (x1 - x0)/(u - v);
        x0 = x1;
        v = u;
        x1 -= u*s;
        u = f(x1);
        if (fabs(x1 - x0) < delta || fabs(u) < eps)
            break;
    }
    return x1;
}
